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LARGE SAMPLE THEORY OF INTRINSIC AND EXTRINSIC 
SAMPLE MEANS ON MANIFOLDS II 

By Rabi Bhattacharya 1 and Vic Patrangenaru 2 
University of Arizona and Texas Tech University 

This article develops nonparametric inference procedures for esti- 
mation and testing problems for means on manifolds. A central limit 
theorem for Frechet sample means is derived leading to an asymptotic 
distribution theory of intrinsic sample means on Riemannian mani- 
folds. Central limit theorems are also obtained for extrinsic sample 
means w.r.t. an arbitrary embedding of a differentiable manifold in 
a Euclidean space. Bootstrap methods particularly suitable for these 
problems are presented. Applications are given to distributions on 
the sphere S d (directional spaces), real projective space RP^ -1 (ax- 
ial spaces), complex projective space CP k ~ 2 (planar shape spaces) 
w.r.t. Veronese-Whitney embeddings and a three-dimensional shape 
space E|. 

1. Introduction. Statistical inference for distributions on manifolds is 
now a broad discipline with wide-ranging applications. Its study has gained 
momentum in recent years, especially due to applications in biosciences and 
medicine, and in image analysis. Including in the substantial body of litera- 
ture in this field are the books by Bookstein [10], Dryden and Mardia [15], 
Kendall, Barden, Carne and Le [33], Mardia and Jupp [41], Small [49] and 
Watson [52]. While much of this literature focuses on parametric or semi- 
parametric models, the present article aims at providing a general framework 
for nonparametric inference for location. This is a continuation of our ear- 
lier work [7, 8] where some general properties of extrinsic and intrinsic mean 
sets on general manifolds were derived, and the problem of consistency of the 
corresponding sample indices was explored. The main focus of the present 
article is the derivation of asymptotic distributions of intrinsic and extrinsic 
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sample means and confidence regions based on them. We provide classical 
CLT-based confidence regions and tests based on them, as well as those 
based on Efron's bootstrap [17]. 

Measures of location and dispersion for distributions on a manifold M 
were studied in [7, 8] as Frechet parameters associated with two types of 
distances on M. If j : M — > M. k is an embedding, the Euclidean distance 
restricted to j(M) yields the extrinsic mean set and the extrinsic total vari- 
ance. On the other hand, a Riemannian distance on M yields the intrinsic 
mean set and intrinsic total variance. 

Recall that the Frechet mean of a probability measure Q on a complete 
metric space (M,p) is the minimizer of the function F(x) = J p 2 (x,y)Q(dy), 
when such a minimizer exists and is unique [21]. In general the set of min- 
imizers of F is called the Frechet mean set. The intrinsic mean fJ,i(Q) is 
the Frechet mean of a probability measure Q on a complete d-dimensional 
Riemannian manifold M endowed with the geodesic distance d g determined 
by the Riemannian structure g on M. It is known that if Q is sufficiently 
concentrated, then /Uj(Q) exists [see Theorem 2.2(a)]. The extrinsic mean 
^e{Q) = Pj,E{Q) of a probability measure Q on a manifold M w.r.t. an em- 
bedding j : M — > M. k is the Frechet mean associated with the restriction to 
j(M) of the Euclidean distance in M. k . In [8] it was shown that the extrinsic 
mean of Q exists if the ordinary mean of j(Q) is a nonfocal point of j(M), 
that is, if there is a unique point xq on j(M) having the smallest distance 
from the mean of j(Q). In this case /J-j t E(Q) = j _1 ( x o)- 

It is easier to compute the intrinsic mean if the Riemannian manifold has 
zero curvature in a neighborhood containing suppQ [45]. In particular this 
is the case for distributions on linear projective shape spaces [42]. If the 
manifold has nonzero curvature around suppQ, it is easier to compute the 
extrinsic sample mean. It may be pointed out that if Q is highly concentrated 
as in our medical examples in [8] and in Section 5, the intrinsic and extrinsic 
means are virtually indistinguishable. 

We now provide a summary of the main results in this article. Section 2 
is devoted to nonparametric inference for the Frechet mean of a probability 
measure Q on a manifold M for which there is a domain U of a chart 
: U — > M. d such that Q(U) = 1. In Theorem 2.1 it is shown that in this case, 
under some rather general assumptions, the image of the Frechet sample 
mean under <f> is asymptotically normally distributed around the image of 
the Frechet mean of Q. This leads to the asymptotic distribution theory 
of the intrinsic sample mean on a Riemannian manifold M (Theorems 2.2, 
2.3). In Corollaries 2.3 and 2.4 bootstrap confidence regions are derived for 
the Frechet mean, with or without a pivot. 

Section 3 is devoted to asymptotics of extrinsic sample means. The ideas 
behind the main result here are essentially due to Hendriks and Landsman 
[27] and Patrangenaru [44]. The two approaches are somewhat different. 
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We present in this article an extension of the latter approach. Extrinsic 
means are commonly used in directional, axial and shape statistics. In the 
particular case of directional data analysis, that is, when M = S^ 1 is the 
unit sphere in M. d , Fisher, Hall, Jing and Wood [19] provided an approach for 
inference using computationally efficient bootstrapping which gets around 
the problem of increased dimensionality associated with the embedding of 
the manifold M in a higher-dimensional Euclidean space. In Corollary 3.2 
confidence regions are derived for the extrinsic mean fJ>j t E(Q)- Nonparametric 
bootstrap methods on abstract manifolds are also derived in this section 
(Theorem 3.2, Proposition 3.2). 

If one assumes that Q has a nonzero absolutely continuous component 
with respect to the volume measure on M, then from some results of Bhat- 
tacharya and Ghosh [6], Babu and Singh [1], Beran [2] and Hall [24, 25], one 
derives bootstrap-based confidence regions for /j>e(Q) with coverage error 
O p (n -2 ) (Theorem 3.4) (also see [5, 9]). One may also use the nonpivotal 
bootstrap to construct confidence regions based on the percentile method of 
Hall [25] for general Q with a coverage error no more than O p (n~ d ^ d+1 ^), 
where d is the dimension of the manifold (see Remark 2.4 and Proposition 
3.2). This is particularly useful in those cases where the asymptotic disper- 
sion matrix is difficult to compute. 

Section 4 applies the preceding theory to (i) real projective spaces IR^" 1 — 
the axial spaces, and (ii) complex projective spaces CP k ~ 2 — the shape spaces. 
Another application to products of real projective spaces (Ri- ,m ) fc_m_1 , or 
the so-called projective shape spaces, will appear in [42]. 

As an application of Corollary 3.3, large sample confidence regions for 
mean axes are described in Corollary 4.2. A similar application to projective 
shape spaces, combining bootstrap methods for directional data from [3], 
appears in [42]. Other applications to axial spaces are given in Theorem 4.3 
and Corollary 4.4, and to planar shape spaces in Theorem 4.5. 

Finally in Section 5 we apply the results of Sections 2 and 4 to construct 
(1) a 95% large sample confidence region for the intrinsic mean location 
of the magnetic South Pole from a directional data set given in [20], (2) 
simultaneous confidence intervals for the affine coordinates of the extrinsic 
sample mean shape in a medical application and (3) a test for the difference 
between three-dimensional mean shapes in a glaucoma detection problem. 

2. A central limit theorem for Frechet sample means and bootstrap- 
ping. A (i-dimensional differentiable manifold is a locally compact sepa- 
rable Hausdorff space M, together with an atlas Am comprising a family 
of charts (U a ,cp a ) of open sets U a covering M, and for each a a home- 
omorphism cf) a of U a onto an open subset of M. d for which the transition 
maps 4> a ■ (f)^ 1 : (f>p(U a n Up) — ► (f> a (U a D Up) are of class C°°. The sets U a 
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are often called coordinate neighborhoods. One may show that a differen- 
tiable manifold is metrizable. We briefly recall some basic notion associated 
with Riemannian manifolds. For details the reader may refer to any stan- 
dard text on differential geometry (e.g., [13, 26], or [38]). A Riemannian 
metric g on a differentiable manifold M is a C°° symmetric positive defi- 
nite tensor field of type (2,0), that is, a family of inner products g p = (■, -) p 
on the tangent spaces T p M,p € M, which is differentiable w.r.t. p. A Rie- 
mannian manifold M is a connected differentiable manifold endowed with 
a Riemannian metric g. The distance p g induced by g is called the geodesic 

distance. For M, p g (p, q) is the infimum of lengths (x(t), x(t)) l J(t\ dt of 

all C 1 -curves x(t),a <t <b, with x(a) =p,x(b) = q. The minimizer satisfies 
a variational equation whose solution is a geodesic curve. There is a unique 
such geodesic curve t — ► j(t) for any initial point 7(0) =p and initial tangent 
vector 7(0) = v. A classical result of Hopf and Rinow states that (M,p g ) is 
complete as a metric space if and only if (M,g) is geodesically complete [i.e., 
every geodesic curve j(t) is defined for all t E R]. These two equivalent prop- 
erties of completeness are in turn equivalent to a third property: all closed 
bounded subsets of (M,p g ) are compact ([13], pages 146 and 147). 

Given q S M, the exponential map Exp ? :U — » M is defined on an open 
neighborhood U of G T q M by the correspondence v — ► 7«(1), where J v (t) 
is the unique geodesic satisfying 7(0) = (7,7(0) = v, provided j(t) extends at 
least to t = 1. Thus if (M, g) is geodesically complete or, equivalently, (M, p g ) 
is complete as a metric space, then Exp g is defined on all of T q M. In this 
article, unless otherwise specified, all Riemannian manifolds are assumed to 
be complete. 

Note that if 7(0) =p and j(t) is a geodesic, it is generally not true that 
the geodesic distance between p and q = 7(^1), say, is minimized by 7(i), < 
t < ti (consider, e.g., the great circles on the sphere S 2 as geodesies). Let 
*o = to(p) be the supremum of all t± > for which this minimization holds. 
If to < 00, then 7 (to) is the cut point of p along 7. The cut locus C{p) of p 
is the union of all cut points of p along all geodesies 7 starting at p [e.g., 
C(p) = {-p} on S 2 }. 

In this article we deal with both intrinsic and extrinsic means. Hence we 
will often consider a general distance p on a differentiable manifold M, but 
assume that (M,p) is complete as a metric space. We consider only those 
probability measures Q on M for which the Frechet mean pjr = pj^(Q) exists. 
Moreover we assume that there is a chart (U,4>) such that Q{U) = 1, and 
PF^U. 

Remark 2.1. The assumption above on the existence of a chart (U,(p) 
such that Q(U) = 1 is less restrictive than it may seem. If g is a Riemannian 
structure on M and Q is absolutely continuous w.r.t. the volume measure, 
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then, for any given p, the complement U of the cut locus C(p) is the domain 
of definition of such a local coordinate system (the coordinate map being the 
inverse of Exp p , the exponential map at p) (see [38], page 100, for details). 

Example 2.1. For the d-dimensional unit sphere, M = S d = {p £ R d+1 :\\p 
1}, with the Riemannian metric induced by the Euclidean metric on R d+1 , 
the exponential map at a given point p € S d is defined on the tangent space 
T p M and is given by 

(2.1) Exp p (v)=cos(\\v\\)p + sm(\\v\\)\\v\\- 1 v 

If x G S d ,x —p, then there is a unique vector u £ T p M such that x = 
Exp p u, and we will label this vector by u = Log p x. Since T p S d = 
v . p = o}, it follows that 

(2.2) Log p x = (1 — (p ■ x) 2 )" 1 / 2 arccos(p ■ x)(x — (p ■ x)p). 

In particular, for d = 2 we consider the orthobasis e±(p), e2{p) £ T p S 2 , where 
P = (P1,P2,P 3 )* € S 2 \{iV, 5} [iV = (0, 0, 1), 5 = (0,0, -1)]: 

ei(p) = ((pi) 2 + fe) 2 )~ 1/2 (-P2,Pi,0)', 

(2.3) e 2 (p) = (-((Pi) 2 + (P2) 2 r 1/2 PW3, 

-(x 2 + y 2 r 1/2 P2 P3 ,((Pi) 2 + (P2) 2 ) 1/2 ) t . 

The logarithmic coordinates of the point x = (xi,X2-,x^) T are given in this 
case by 

/ 24 ) ^(P) = ei(p) -LogpX, 

^ 2 (p) = e2(p) -LogpX. 

For computations one may use a-b = a l b. 

Now the image measure of Q under (f> has the Frechet mean fi = 4>(pf) 
w.r.t. the distance p^(u,v) := p(4>~ 1 (u) , (ft -1 i v ) )i u i v S <j>(U). Similarly, if X/ 
(i = 1, . . . , n) are i.i.d. with common distribution Q and defined on a prob- 
ability space (fi,.A, P), let /u n] jr be a measurable selection from the Frechet 
mean set (w.r.t. p) of the empirical Q n = ~Y^i=i^Xi- Then p n = 4>(p n ^) 
is a measurable selection from the Frechet mean set (w.r.t. p^) of Qt, = 
— Ya=i$x .) where X, = 0(Xj). Assuming twice continuous differentiability 
of — > {p^) 2 {u,9) 7 write the Euclidean gradient as 

(2.5) *M)=grad (/) 2 (^)= (^^(u^))^ ^(^(ujfl))^. 
Now p is the point of minimum of 

(2.6) F+(0):= J(p*) 2 (u,9)Q^du) 



G 



R. BHATTACHARYA AND V. PATRANGENARU 



and fi n is a local minimum of 

Therefore, one has the Taylor expansion 
1 n 



n 1— ' 

i=i 



1 n 



n d 

+ - E E D r ,^(X t ;^(^ - //') + K (1 < r < d), 



i=l r'=l 

where 



d -i n 



(2.8) R r n = ^^-/l-EP^'^;^)-^^^^)} 



n . 

r'=l i=l 



and 6> n lies on the line segment joining // and \x n (for sufficiently large n). 
We will assume 



r?q x £|*(A^)| 2 <oc, 

1 J £|A.,* r (l i; /x)| 2 <oo (Vr,r'). 

To show that R r n is negligible, write 



u r ' r (x,e) := sup |Z? r '^ r (x; 0) - Z? r /^ r (a;; fj)\ 

and assume 

(2.10) 5 r ' r '(c):=Eu r ' r '(Xi,c)^0 as c j (1 < r, r' < d). 
One may then rewrite (2.7) in vectorial form as 

1 n 

(2.11) = -= *(Xi;») + (A + x/^Mn - H), 
where 

(2-12) A = E((D r ,* r (X l ; f ,))) d rr , =1 

and <5 n — > in probability as n — > oo, if ^ n — ► /i in probability. If, finally, we 
assume A is nonsingular, then (2.11) leads to the equation 

(2.13) v ^(^ n -/x)=A- 1 ^^(A > 4 ; / ^ +6' n , 

where 8' n goes to zero in probability as n — ► oo. We have then arrived at the 
following theorem. 
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Theorem 2.1 (CLT for Frechet sample means). Let Q be a probability 
measure on a differentiable manifold M endowed with a metric p such that 
every closed and bounded set of (M,p) is compact. Assume (i) the Frechet 
mean p^ exists, (ii) there exists a coordinate neighborhood (U, eft) such that 
Q(U) = 1, (iii) the map 6 — > (p^) 2 (6,u) is twice continuously differentiable on 
4>{U), (iv) the integrability conditions (2.9) hold as well as the relation (2.10) 
and (v) A, defined by (2.12), is nonsingular. Then (a) every measurable 
selection p n from the (sample) Frechet mean set of = -J2i=i°~x. * s a 

consistent estimator of p, and (b) yjn{p n — p) — > Af(0, A~ 1 C(A')~ 1 ), where 
C is the covariance matrix of ^{Xi]p). 

Proof. Part (a) follows from Theorem 2.3 in [8]. The proof of part (b) 
is as outlined above, and it may also be derived from standard proofs of the 
CLT for M-estimators (see, e.g., [29], pages 132-134). □ 

As an immediate corollary one obtains: 

Corollary 2.1. Let (M,g) be a Riemannian manifold and let p = p g be 
the geodesic distance. Let Q be a probability measure on M whose support 
is compact and is contained in a coordinate neighborhood (U,(j>). Assume 
that (i) the intrinsic mean pj = pj? exists, (ii) the map 9 — > (p^) 2 (9,u) is 
twice continuously differentiable on 4>{U) for each uE 4>(U) and A, defined 
by (2.12), is nonsingular. Then the conclusions of Theorem 2.1 hold for the 
intrinsic sample mean p n j = p n ^ of Q n = - Ya=i $Xh with p = <fi{pi). 

We now prove one of the main results of this section. 

Theorem 2.2 (CLT for intrinsic sample means). Let (M,g) be a Rie- 
mannian manifold and let p = p g be the geodesic distance. Let Q be a prob- 
ability measure on M whose support is contained in a closed geodesic ball 
B r = B r (xo) with center xo and radius r which is disjoint from the cut locus 
C{xq). Assume r < where K 2 is the supremum of sectional curvatures 
in B r if this supremum is positive, or zero if this supremum is nonposi- 
tive. Then (a) the intrinsic mean pj (of Q) exists, and (b) the conclusion 
of Theorem 2.1 holds for the image p n = 4>{p n ,l) of the intrinsic sample 
mean p n j of Q n = - Ya=i °~x , under the inverse 4> of the exponential map, 
<f>=(ExpJ-\ 

Proof, (a) It is known that under the given assumptions, there is a 
local minimum pj, say, of the Frechet function F which belongs to B r and 
that this minimum is also the unique minimum in Z?2r [30, 34, 40]. We now 
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show that pi is actually the unique global minimum of F. Let p E (B2 r ) c - 
Then p{p,x) > r, V x E B r . Hence 

(2.14) F(p)= f_ p 2 (p,x)Q{dx)> f_ r 2 Q(dx) = r 2 . 
On the other hand, 

(2.15) F(mi) < F(xo) = [_ p 2 (x ,x)Q(dx) < r 2 , 

proving F(p) > F(/j,j). 

(b) In view of Corollary 2.1, we only need to show that the Hessian matrix 
A = A(/i) of F o at fx := 4>{pi) is nonsingular, where (f> = Exp" 1 . Now 
according to [30], Theorem 1.2, for every geodesic curve ^(t) in B r ,t E (c,d) 
for some c < 0, d > 0, 

(2.16) ^F( 7 (i))>0 (c<t<d). 

Let ip = Exp M/ denote the exponential map at /i/, and let 7(t) be the unique 
geodesic with 7(0) = pi and 7 (0) = v, so that 7 (f) = tp(tv). Here we identify 
the tangent space T^M with R d . Applying (2.16) to this geodesic (at t = 0), 
and writing G = F o one has 



(2.17) ^H^tv)) 



= J2v i v j (D i D j G){0)>0 (Vv^O), 

t=o 



that is, the Hessian of G is positive definite at E If xq = pj, this 
completes the proof of (b). 

Next let xq 7^ /U/. Now Fo(j)~ l = Go (ip -1 o cj)" 1 ) on a domain that includes 
p = 4>(fJ>i) = (Exp a , )~ 1 (/i/). Write f^" 1 o 0" 1 = /. Then in a neighborhood 
of fi, 

9 2 {Gof) l _ ap a/j" 

+EPiG)c/(«))^(«)- 

The second sum in (2.18) vanishes at u = p,, since (DjG)(f(p)) = (DjG)(0) = 
as f(p,) = ip" 1 4>~ l (p) = ^~ l {pi) = is a local minimum of G. Also / is a 
diffeomorphism in a neighborhood of p. Hence, writing A rr r(p) as the (r, r') 
element of A(/i), 

KM = dy °/ ^ (u)=Y(D j D i ,G)(0)^(ji)^(jj). 
r,r\v; du r du r x ' z - 1 , 3 n 'du r ^'du r K ' 

3,3 

This shows, along with (2.17), that A = A(/z) is positive definite. □ 
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Remark 2.2. If the supremum of the sectional curvatures (of a complete 
manifold M) is nonpositive, and the support of Q is contained in B r , then 
the hypotheses of Theorem 2.2 are satisfied, and the conclusions (a), (b) 
hold. One may apply this even with r = oo. 

Remark 2.3. The assumptions in Theorem 2.2 on the support of Q for 
the existence of \i\ are too restrictive for general applications. But without 
additional structures they cannot be entirely dispensed with, as is easily 
shown by letting Q be the uniform distribution on the equator of S 2 . For the 
complex projective space even, necessary and sufficient conditions 

for the existence of the intrinsic mean fij of an absolutely continuous (w.r.t. 
the volume measure) Q with radially symmetric density are given in [33, 39]. 

It may be pointed out that it is the assumption of some symmetry, that 
is, the invariance of Q under a group of isometries, that often causes the in- 
trinsic mean set to contain more than one element (see, e.g., [8], Proposition 
2.2). The next result is, therefore, expected to be more generally applicable 
than Theorem 2.2. 



Theorem 2.3 (CLT for intrinsic sample means). Let Q be absolutely 
continuous w.r.t. the volume measure on a Riemannian manifold (M,g). 
Assume that (i) fij exists, (ii) the integrability conditions (2.9) hold, (iii) the 
Hessian matrix A of F o cp" 1 at fi = (p(fJ-i) is nonsingular and (iv) the co- 

variance matrix C of '\P \Xf, fi) is nonsingular. Then y / n(/x n — fj) — >Af(0,T), 
where F = A _1 C(A*) _1 . 

This theorem follows from Theorem 2.1 and Remark 2.1. 

In order to obtain a confidence region for jijr using the CLT in Theo- 
rem 2.1 in the traditional manner, one needs to estimate the covariance 
matrix F = A C(A*) _1 . For this one may use proper estimates of A and C, 
namely, 

1 n 

(2 19) K6) -E( Grad *X X -/^)> C = CovQt 

f ■-A- 1 C(A t )~ 1 , f-^A'C- 1 !. 

The following corollary is now immediate. Let Xdi-a denote the (1 — a)th 
quantile of the chi-square distribution with d degrees of freedom. 

Corollary 2.2. Under the hypothesis of Theorem 2.1, if C is nonsin- 
gular, a confidence region for [i? of asymptotic level 1 — a is given by U nt0l := 
^ _1 (-D n ,a), where D n a = {v G <j){U) : n(fj, n - t>)*f _1 (/i n - v) < xli-al- 
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Example 2.2. In the case of the sphere S 2 from Example 2.1, it fol- 
lows that if we consider an arbitrary data point u = {u l ,u 2 ), and a second 
point 6 = Log„A = (6 l ,9 2 ), and evaluate the matrix of second-order partial 
derivatives w.r.t. 8 l ,6 2 of 

(2.20) G(u, 6) = arccos 2 ( cos llull + sm ^ (u 1 6 1 + u 2 d 2 ) - -\\9\\ 2 cos llul 

V IMI 2 

then 

d 2 G . , 2u r u s ( llull \ 2$„||u|| 

( 2 - 21 ) *™ («; o) = w 1 - 73 + 



d6 r 80 s ' IMI 2 V tan||u||/ tan||u||' 

where <5 rs is the Kronecker symbol and ||u|| 2 = (u 1 ) 2 + (u 2 ) 2 . The matrix 
A = (Xrr')r,r'=i,2 nas the entries 

1 n d 2 G 

2 = 1 

Assume C is the sample covariance matrix of Uj, j = 1, . . . ,n; a large sample 
confidence region for the intrinsic mean is given by Corollary 2.2 with fi n = 0. 



We now turn to the problem of bootstrapping a confidence region for 
fip. Let X* n be i.i.d. with common distribution Q n (conditionally, given 
{Xi :l<i< n}). Write X* n = <p(X* ), l<i<n, and let /x* be a measurable 
selection from the Frechet mean set of Q*'^ := ^Ya=i^x* ■ Let £* a be a 
subset of 4>(U), such that P*(/x* — fi n £ E^ a ) — > 1 — a in probability, where 
P* denotes the probability under Q n . 



Corollary 2.3. In addition to the hypothesis of Theorem 2.1, assume 
C is nonsingular. Then (f)" 1 ({(fi n — E* a ) C\ (p(U)}) is a confidence region for 
HP of asymptotic level (1 — a). 

Proof. One may write (2.7) and (2.8) with // and fx n replaced by 
/i n and /i*, respectively, also replacing Xi by X* in (2.8). To show that a new 
version of (2.11) holds with similar replacements (also replacing A by A), 
with a 5^ (in place of 5 n ) going to zero in probability, one may apply Cheby- 
shev's inequality with a first-order absolute moment under Q n , proving 
that A* - A goes to zero in probability. Here A* = \ V^^Grad*)^;^*). 
One then arrives at the desired version of (2.7), replacing fj, n , /j,,A,Xi by 
//*,// n , A, X*, respectively, and with the remainder (corresponding to 8 n ) 
going to zero in probability. □ 
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Remark 2.4. In Corollary 2.3, we have considered the so-called per- 
centile bootstrap of Hall [25] (also see [17]), which does not require the com- 
putation of the standard error A. For this as well as for the CLT-based 
confidence region given by Corollary 2.2, one can show that the coverage 
error is no more than O p (n~ d ^ d+1 ^) or 0(n~ d ^ d+1 ^), as the case may be 
[4]. One may also use the bootstrap distribution of the pivotal statistic 
n(/j, n — /j,) T r~ 1 ([i n — fi) to find c* a such that 

(2.23) P*(n(fi* n - finff *~H< ~ Mn) < <J - 1 - «, 
to find the confidence region 

(2.24) D* >a = {ve 4>{U) : n(/x n - vft~ l (fx n - v) < < >Q }. 

In particular, if Q has a nonzero absolutely continuous component w.r.t. the 
volume measure on M, then so does w.r.t. the Lebesgue measure on 4>(U) 
(see [13], page 44). Then assuming (a) c* a is such that the P*-probability 
in (2.23) equals 1 — a + O p (n~ 2 ) and (b) some additional smoothness and 
integrability conditions of the third derivatives of \E r , one can show that the 
coverage error [i.e., the difference between 1 — a and P(/j, € -D* Q )] is O p (n~ 2 ) 
(see [5, 6, 12, 24, 25]). It follows that the coverage error of the confidence 
region <^> _1 (Z)* a n <fi(U)) for fijr is also 0(n~ 2 ). We state one such result 
precisely. 

Corollary 2.4 (Bootstrapping the intrinsic sample mean). Suppose 
the hypothesis of Theorem 2.3 holds. Then 

sup|P*(ra« - /J, n ) T f - /in) < r) 

- P(n(n n - nff- 1 ^ -fj,)<r)\ = O p (n~ 2 ), 
and the coverage error of the pivotal bootstrap confidence region is = O p (n~ 2 ). 

Remark 2.5. The assumption of absolute continuity of Q in Theorem 
2.3 is reasonable for most applications. Indeed this is assumed in most para- 
metric models in directional and shape analysis (see, e.g., [15, 52]). 

Remark 2.6. The results of this section may be extended to the two- 
sample problem, or to paired samples, in a fairly straightforward manner. 
For example, in the case of paired observations (Xi,Yi),i = 1, . . . ,n, let X{ 
have (marginal) distribution Q, and intrinsic mean /i/, and let Q2 and vj 
be the corresponding quantities for Y\ . Let <f> = Exp^ 1 for some xq , and let 
fi, v and n n ,v n be the images under <j) of the intrinsic population and sample 
means. Then one arrives at the following [see (2.13)]: 



(2.25) 



Vn(fi n - (j,)- ^fn{y n - v) -> JV(0, V), 
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where T is the covariance matrix of Aj~ 1 ^(Xj; y) — A^" 1 \I'(1^; v). Here Aj is 
the Hessian matrix of F o cf)~ l for Qi (i = 1,2). Assume T is nonsingular. 
Then a CLT-based confidence region for 7 := f/, — v is given in terms of 
In ■= ~ v n by {v £ R d : n(7 n - v)f _1 (7,, - «) < Xd,i_ a }- Alternatively, one 
may use a bootstrap estimate of the distribution of \fn{^ n — 7) to derive a 
confidence region. 

In Section 5 we consider two applications of results in this section (and 
one application of the results in Sections 3 and 4). Application 1 deals with 
the data from a paleomagnetic study of the possible migration of the Earth's 
magnetic poles over geological time scales. Here M = S 2 and the geodesic 
distance between two points is the arclength between them measured on the 
great circle passing through them. 

Application 3 analyzes some recent three-dimensional image data on the 
effect of a (temporary) glaucoma-inducing treatment in 12 Rhesus monkeys. 
On each animal k = 4 carefully chosen landmarks are measured on each 
eye — the normal eye and the treated eye. For each observation (a set of 
four points in M 3 ) the effects of translation, rotation and size are removed 
to obtain a sample of 12 points on the five-dimensional shape orbifold S|. 
We use the so-called three-dimensional Bookstein coordinates to label these 
points (see [15], pages 78-80). In order to apply Theorem 2.3 (i.e., its analog 
indicated above), a somewhat flat Riemannian structure is chosen so that 
the necessary assumptions can be verified. 

3. The CLT for extrinsic sample means and confidence regions for the 
extrinsic mean. From Theorem 2.1 one may derive a CLT for extrinsic 
sample means similar to Corollary 2.1. In this section, however, we use an- 
other approach which, for extrinsic means, is simpler to apply and generally 
less restrictive. 

Recall that the extrinsic mean Hj^iQ) of a nonfocal probability measure 
Q on a manifold M w.r.t. an embedding j : M — > W k , when it exists, is given 
by Hj,E{Q) = 3 {Pj{^))t where fj, is the mean of j(Q) and Pj is the projec- 
tion on j(M) (see [8], Proposition 3.1, e.g.). Often the extrinsic mean will 
be denoted by he{Q), or simply he, when j and Q are fixed in a particular 
context. To ensure the existence of the extrinsic mean set, in this section we 
will assume that j(M) is closed in W k . 

Assume (X\, . . . ,X n ) are i.i.d. M- valued random objects whose common 
probability distribution is Q, and let Xe '■= He{Qti) be the extrinsic sample 
mean. Here Q n = - Y7j=i $Xj is the empirical distribution. 

A CLT for the extrinsic sample mean on a submanifold M of M. k (with j 
the inclusion map) was derived by Hendriks and Landsman [27] and, inde- 
pendently, by Patrangenaru [44] by different methods. Differentiable man- 
ifolds that are not a priori submanifolds of M. k arise in new areas of data 
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analysis such as in shape analysis, in high-level image analysis, or in signal 
and image processing (see, e.g., [15, 16, 22, 31, 32, 33, 42, 51]). These man- 
ifolds, known under the names of shape spaces and projective shape spaces, 
are quotient spaces of submanifolds of M k (spaces of orbits of actions of Lie 
groups), rather than submanifolds of IR fe . Our approach is a generalization of 
the adapted frame method of Patrangenaru [44] to closed embeddings in 
This method leads to an appropriate dimension reduction in the CLT and, 
thereby, reduces computational intensity. This method extends the results 
of Fisher et al. [19] who considered the case M = S d . We expect that with 
some effort the results of Hendriks and Landsman [27] may be modified to 
yield the same result. 

Assume j is an embedding of a d-dimensional manifold M such that 
j(M) is closed in M. k , and Q is a j-nonfocal probability measure on M 
such that j(Q) has finite moments of order 2 (or of sufficiently high order 
as needed). Let [i and £ be, respectively, the mean and covariance ma- 
trix of j(Q) regarded as a probability measure on R fc . Let T be the set of 
focal points of j(M), and let Pj : .P c — > j(M) be the projection on j(M). 
Pj is differentiable at fj, and has the differentiability class of j{M) around 
any nonfocal point. In order to evaluate the differential d^Pj we consider 
a special orthonormal frame field that will ease the computations. Assume 
p — > (f\(x), . . . , fd{x)) is a local frame field on an open subset of M such that, 
for each x G M, (dj(f\(x)), . . . ,dj(fd{x))) are orthonormal vectors in M. k . A 
local frame field (ei (p) , e 2 (p) , . . . , (p) ) defined on an open neighborhood 
U C M fc is adapted to the embedding j if it is an orthonormal frame field 
and V.x £ j~ 1 (C7), (e r (j(x)) = d p j(f r (x)), r = l,...,d. Let ei,e 2 ,...,e fc be 
the canonical basis of M fc and assume (e\(p), e2(p), . . . , e^(p)) is an adapted 
frame field around Pj(n) = j{he)- Then d^Pj(eb) G T P .^j(M) is a linear 
combination of ei(Pj (//)), e2(Pj(/x) ),..., ed(Pj{p))\ 

(3.1) d^P ] {e b ) = Y J {d^P ] {e b ))-e a {P 3 {^))e a (P ] { l x)). 

By the delta method, n 1 / 2 (P j (j(X)) - converges weakly to N(0, £ ), 

where = ± E? =1 J« and 



(3-2) % 



=i 

[^d M P j (e 6 )- ea (P : ,(/x))e a (P i ( M )) 



a=l 
X 



V 



b=l....,k 
t 



6=1, 



Here £is the covariance matrix of j(X\) w.r.t. the canonical basis e±, e2, • • • , &k- 
The asymptotic distribution N(0, £ ) is degenerate and can be regarded as 
a distribution on Tp ^j(M), since the range of d^Pj is T P ^j(M). Note 
that 

d^Pj(e b ) ■ e a (Pj(n)) = for a = d+ 1, . . . , k. 
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Remark 3.1. An asymptotic distribution of the extrinsic sample mean 
can be obtained as a particular case of Theorem 2.1. The covariance matrix 
in that theorem depends both on the way the manifold is embedded and 
on the chart used. We provide below an alternative CLT, which applies to 
arbitrary embeddings, leads to pivots and is independent of the chart used. 

The tangential component tan(f) of v E M. k w.r.t. the basis e a (Pj(^)) G 
T Pj{n)i( M ), a = 1> • • • > d, is given by 

(3.3) tan(u) = ( ei (P»)V e d (P 3 (^)) T v) 1 . 

Then the random vector (d tlE j)~ 1 {taii(Pj(j (X)) — Pj(fi))) = Ylt,=i ~Xjfa nas 
the following covariance matrix w.r.t. the basis fi(fJ-E), ■ ■ ■ , fd{^E)'- 

Zj,E = e (-Pj(^))*S M e 6 (P j (/i)) 1<a6 < d 

(3.4) = \j2 d » P i( e b) ' ea ( P i^))\ a=l td V 



d^Pj{e b ) ■ e a (Pj(n)) 



t 



0=1, 



Definition 3.1. The matrix given by (3.4) is the extrinsic co- 
variance matrix of the j-nonfocal distribution Q (of X\) w.r.t. the basis 
fl(jiE),---,fd(V>E)- 

When j is fixed in a specific context, the subscript j in will be 

omitted. If, in addition, rank t = d, £j,.e is invertible and we define the 
j- standardized mean vector 

(3.5) Z j , n =:n i /%, E -V\X 1 j ,...,X d j f. 

Proposition 3.1. Assume {X r } r= i n is a random sample from the 
j-nonfocal distribution j(Q), and let fj, = E(j(X\)) and assume the extrin- 
sic covariance matrix E^e of Q is finite. Let (ei(p), e2(p), ■ ■ • , e&(p)) be an 
orthonormal frame field adapted to j. Then (a) the extrinsic sample mean 
Xe has asymptotically a normal distribution in the tangent space to M at 
He{Q) with mean and covariance matrix n~ 1 T,j t E, and (b) ifSj^E is non- 
singular, the j -standardized mean vector Z^ n given in (3.5) converges weakly 
to N(0,I d ). " 

As a particular case of Proposition 3.1, when j is the inclusion map of a 
submanifold of M. k , we get the following result for nonfocal distributions on 
an arbitrary closed submanifold M of M. k : 

Corollary 3.1. Assume M C R k is a closed submanifold ofR k . Let 
{X r } r= i ... n be a random sample from the nonfocal distribution Q on M , 
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and let [X = E(X\) and assume the covariance matrix £ of j(Q) is finite. 
Let (e±(p), e2(p), . . . , e&(p)) be an orthonormal frame field adapted to M . Let 
T*e '■= Tjj,Ei where j : M — > M k is the inclusion map. Then (a) n 1 / 2 tan.(J(X e) — 
Hhe)) converges weakly to N(0,T,e), and (b) if$ induces a nonsingular bi- 
linear form on Tj^ E ^j{M), then ||Zj n || 2 converges weakly to the chi-square 
distribution x\- 

Example 3.1. In the case of a hypersphere in j(x) = x and P, = 
Pm- We evaluate the statistic ||Zj n || 2 = njlE^ -1 / 2 tan(Pw(X) — Pm(^))\\ 2 ■ 
The projection map is Pm{x) =x/||x||. Pm has the following property: if 
v = cx, then d x PM(v) = 0; on the other hand, if the restriction of d x Pyi to 
the orthocomplement of Kx is a conformal map, that is, if v ■ x = 0, then 
dxPuiv) = . In particular, if we select the coordinate system such 

that x = ||x||efc, then one may take e a (PM(x)) = e a , and we get 

d x P M (e b ) ■e a (P M {x)) = INI" 1 ^ Va,6=l,...,fe- l,d x P M (e k ) = 0. 

Since efc(PM (/■*)) points in the direction of /i, c^-Pm^) • fJ- = 0, V6 = 1, . . . , k — 
1, and we get 

(3.6) E E = ||^ir 2 ^([X • e (/VllM||)]a=i,...,fc-i[X • e ( M /||/i||)]* = i,..., fc -i) 
which is the matrix G in formula (A.l) in [19]. 



Remark 3.2. The CLT for extrinsic sample means as stated in Propo- 
sition 3.1 or Corollary 3.1 cannot be used to construct confidence regions 
for extrinsic means, since the population extrinsic covariance matrix is un- 
known. In order to find a consistent estimator of Ej,E, note that j(X) is 
a consistent estimator of fx, djr^r Pj converges in probability to d^Pj, and 

e a (Pj(J(X))) converges in probability to e a {Pj{fx)) and, further, 

Sj,n = n~ l J2ti(Xr) -P))(j'(lr) " P)f 

is a consistent estimator of £L It follows that 



(3.7) 



J2 d m Pj{e b ) ■ e a (P 3 (j(X)))e a (P 3 (j(X))) 



la=l 



S 



dj^Pjieb) ■ e a {P 3 {j{X)))e a {P 3 {j{X))) 



a=l 



is a consistent estimator of S^, and tan p .(jfxj) 



v is a consistent estimator of 



tan(u). 
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If we take the components of the bilinear form associated with the ma- 
trix (3.7) w.r.t. ei (Pj0(Tj)), e 2 (P ? (j(X))), . . . , e d (Pj0(Tj)), we get a con- 
sistent estimator of ^j,e given by 



G(j,X) = [[J2<%n P j(ei>)-e*( p MX))) 

x [[E^iW-e a (Pi(P))) 



a=l,...,a 
0=1, 



and obtain the following results. 



Theorem 3.1. Assume j:M^R k is a closed embedding of M in R k . 
Let {X r } r 

=l....,n be a random sample from the j -nonfocal distribution Q, and 
let fi = E(j(Xi)) and assume j{X\) has finite second- order moments and the 
extrinsic covariance matrix of X\ is nonsingular. Let (ei(p), e2(p), . . . , &k{p)) 
be an orthonormal frame field adapted to j . IfG(j,X) is given by (3.8), then 
for n large enough G(j,X) is nonsingular (with probability converging to 1) 
and (a) the statistic 

(3.9) n 1 ' 2 GV,X)-VHan{P j (j(J0)-Pi(t*)) 
converges weakly to N(0,ld), so that 

(3.10) n\\G{j,X)-^HMP3{m)- p mt 
converges weakly to \% an d (b) the statistic 

(3.11) n^GU^y^Han^-^^UiX)) - P») 
converges weakly to N(0,ld), so that 

(3.12) n||G0\X)-V2 tai i p .^(P i (i(X))-P f ( / i))|| 2 
converges weakly to x\- 

Corollary 3.2. Under the hypothesis of Theorem 3.1, a confidence re- 
gion for he of asymptotic level 1 — a is given by (a) C n)Q := j~ 1 (U rita ), where 
U n , a = {fi £ i(M):n||G(i,X)- 1 /2 tan (P i (J(Z)) _ p^f < r by 



(b) D^ a :=r l {V n , a ), whereV n>a = {»ej(M):n\\G(j,X)-V 2 xta^ p , Gm) (P j (j(X))- 

PAp))\\ 2 <xh- a }. 

Theorem 3.1 and Corollary 3.2 involve pivotal statistics. The advantages 
of using pivotal statistics in bootstrapping for confidence regions are well 
known (see, e.g., [1, 2, 5, 9, 24, 25]). 

At this point we recall the steps that one takes to obtain a bootstrapped 
statistic from a pivotal statistic. If {X r } r —i ... n is a random sample from 
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the unknown distribution Q, and {^*} r =i,...,n is a random sample from the 
empirical Q n , conditionally given {X r } r= \ n , then the statistic 

T{X,Q) = n||G(i,X)- 1 /2 t an(P ;7 (i(X)) - P^f 

given in Theorem 3.1(a) has the bootstrap analog 

T(X\Q n )=n\\G(j,X*)- l lH a n P]{1W)) {P 3 ^ 

Here G(j, X* ) is obtained from G(j, X) by substituting X\ , . . . , X* for X±, . . . , X, 
and T(X*,Q n ) is obtained from T(X,Q) by substituting X*, . . . ,X* for 
X!,...,X n , j(Xj) for /i and G(j,X*) for 

The same procedure can be used for the vector-valued statistic 

V(X,Q)=n 1 / 2 G(j,Xy 1 / 2 tMPj6(X))-P J (ri), 
and as a result we get the bootstrapped statistic 

F*(X^Q n )=nV^x*)-V2 tan ^ ( _ )( p. ( ^ ) _p.(^ )) . 

For the rest of this section, we will assume that j(Q), when viewed as 
a measure on the ambient space has finite moments of sufficiently high 
order. If M is compact, then this is automatic. In the noncompact case finite- 
ness of moments of order 12, along with an assumption of a nonzero abso- 
lutely continuous component, is sufficient to ensure an Edgeworth expansion 
up to order 0(n~ 2 ) of the pivotal statistic V(X,Q) (see [5, 6, 12, 19, 24]). 
We then obtain the following results: 

Theorem 3.2. Let{X r } r 

=i,...,n frc a random sample from the j -nonfocal 
distribution Q which has a nonzero absolutely continuous component w.r.t. 
the volume measure on M induced by j . Let fi = E{j{X\)) and assume the 
covariance matrix £ of j (X\ ) is defined and the extrinsic covariance matrix 
is nonsingular and let (ei(p), e2(p), . . . ,ek(p)) be an orthonormal frame 
field adapted to j. Then the distribution function of 

n\\G(j,X)-^tMP 3 (W)) - P^))\\ 2 
can be approximated by the bootstrap distribution function of 

n\\G(j,X*)- 1/2 tan p . ^ ( Pj - Pj (j(X)) ) 1 1 2 

with a coverage error O p (n~ 2 ). 

One may also use nonpivotal bootstrap confidence regions, especially 
when G(j,X) is difficult to compute. The result in this case is the following 
(see [4]). 



18 R. BHATTACHARYA AND V. PATRANGENARU 

Proposition 3.2. Under the hypothesis of Proposition 3.1, the distribu- 
tion function of n\\taxi(Pj(j(X)) — Pj(fi))\\ 2 can be approximated uniformly 
by the bootstrap distribution of 

n||tan p . ( ^ (P^jp^) - P,(j(X)))|| 2 
to provide a confidence region for he with a coverage error no more than 

Remark 3.3. Note that Corollary 3.2(b) provides a computationally 
simpler scheme than Corollary 3.2(a) for large sample confidence regions; but 
for bootstrap confidence regions Theorem 3.2, which is the bootstrap analog 
of Corollary 3.2(a), yields a simpler method. The corresponding 100(1 — a)% 
confidence region is C* a :=j~ 1 (U* a ) with U. * a given by 

Ka = i» G 3(M) :n\\G{j,X)- 1 / 2 tan(P,(j(X)) - P^))f < c^J, 

(3.13) 

where c\_ a is the upper 100(1 — a)% point of the values 

(3.14) «||G(i ) X*)- 1 / 2 tan p . ra (P^j(^)-P i (j(X)))|| 2 

among the bootstrap resamples. One could also use the bootstrap analog 
of the confidence region given in Corollary 3.2(b) for which the confidence 
region is D* a :=i _1 (K* a ) with V* a given by 

K,a = {M e 3( M ) '■ 

^ n\\G{j,Xy 1 / 2 toi Pim) {Pti(X)) - P»)f < dl- a }, 

where d\_ a is the upper 100(1 — a)% point of the values 
(3.16) nllG^X^-V^a^ —^p^^) - Pj(j(X)))\\ 2 

among the bootstrap resamples. The region given by (3.13)-(3.14) has cov- 
erage error O p (n~ 2 ). 

4. Asymptotic distributions of sample mean axes, Procrustes mean shapes 
and extrinsic mean planar projective shapes. In this section we focus on 
the asymptotic distribution of sample means in axial data analysis and 
in planar shape data analysis. The axial space is the (N — l)-dimensional 
real projective space M = RP N ~ l which can be identified with the sphere 
S N ~ l ={x£ K^HI^II 2 = 1} with antipodal points identified (see, e.g., [41]). 
If [x] = {x, —x} G RP^" 1 , \\x\\ = 1, the tangent space at [x] can be described 
as 

(4.1) T^jMP^- 1 = {([x],v),v eR N \v t x = 0}. 
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We consider here the general situation when the distribution on MP 
may not be concentrated. Note that for N odd, RP^" 1 cannot be embedded 
in R , since for any embedding of M.P N ~ l in R fc with N odd, the first Stiefel- 
Whitney class of the normal bundle is not zero ([43], page 51). 

The Veronese-Whitney embedding is defined for arbitrary N by the for- 
mula 

(4.2) j([x])=xx\ \\x\\ = l. 

The embedding j maps M.P N ~ l into a (^N(N + 1) — l)-dimensional Eu- 
clidean hypersphere in the space S(N, R) of real N x N symmetric matrices, 
where the Euclidean distance do between two symmetric matrices is 

d (A,B) = Tii(A-B) 2 ). 

This embedding, which was already used by Watson [52], is preferred over 
other embeddings in Euclidean spaces because it is equivariant (see [35]). 
This means that the special orthogonal group SO(iV) of orthogonal matrices 
with determinant +1 acts as a group of isometries on IP^" 1 with the metric 
of constant positive curvature; and it also acts on the left on S+(N,M), 
the set of nonnegative definite symmetric matrices with real coefficients, by 
T-A = TAT 1 . Also, j(T ■ [x]) = T ■ j([x}), VT e SO(iV), V [x] G RP N ~ 1 . 

Note that j(M.P N ~ 1 ) is the set of all nonnegative definite matrices in 
S(N, R) of rank 1 and trace 1. The following result appears in [8]. 

PROPOSITION 4.1. (a) The set T of the focal points of j(RP N ~ 1 ) in 
<S_l_(iV, R) is the set of matrices in S+(N, R) whose largest eigenvalues are 
of multiplicity at least 2. (b) The projection P, : S+(N, R)\T -> ^RP^" 1 ) 
assigns to each nonnegative definite symmetric matrix A with a highest 
eigenvalue of multiplicity 1, the matrix j([m\), where m(||m|| = 1) is an 
eigenvector of A corresponding to its largest eigenvalue. 

The following result of Prentice [46] is also needed in the sequel. 

Proposition 4.2 ([46]). Assume [X r ], \\X r \\ = 1, r = 1, ... ,n, is a ran- 
dom sample from a j-nonfocal probability measure Q on RP^ -1 . Then the 
j-extrinsic sample covariance matrix G(j,X) is given by 

G{j,X) ab = n~ 1 (r] N -ri a y l (r] N - %) _1 
( 4 - 3 ) xJ2(m a -X r )(m b -X r )(m-X r ) 2 , 

r 

where r] a ,a = l,...,N, are eigenvalues of K := n _1 J2r=i ^rX^ in increas- 
ing order and m a ,a = 1, . . . , N, are corresponding linearly independent unit 
eigenvectors. 
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Here we give a proof of (4.3) based on the equivariance of j to prepare the 
reader for a similar but more complicated formula of the analogous estimator 
given later for CP k ~ 2 . 

Since the map j is equivariant, w.l.o.g. one may assume that j(Xe) = 
Pj(j(X)) is a diagonal matrix, Xe = [m^] = [ejv] and the other unit eigen- 
vectors of j(X) = D are m a = e a , Va = 1, . . . ,N— 1. We evaluate dr>Py Based 
on this description of T^jRP^ -1 , one can select in Tp.^j{M.P N ~ l ) the or- 
thonormal frame e a (Pj(D)) = du N -\j(e a ). Note that S(N, R) has the orthoba- 
sis F%, b <a, where, for a <b, the matrix F% has all entries zero except for 
those in the positions (a, 6), (6, a) that are equal to 2 -1 / 2 ; also F£ = j([e a ]). 
A straightforward computation shows that if rj a ,a = 1,...,N, are the eigen- 
values of D in their increasing order, then d D Pj(F^) = 0, \/b < a < N and 
dnPj{F^) = (r/jy — i]a)~ l e a {Pj(D))', from this equation it follows that, if 
j'PO is a diagonal matrix D, then the entry G{j,X) a b is given by 

G{j,X)a b =n- 1 { m -Va)~\vN -VbT^^XW) 2 - 



Taking j(X) to be a diagonal matrix and m a — e a , (4.3) follows. 

Note that hej = [vn], where (u a ),a = 1, . . . , N, are unit eigenvectors of 
E{XX l ) = E(J(Q)) corresponding to eigenvalues in their increasing order. 
Let T([v\) = n\\G{j,X)- l l 2 t?m{P j (j{X)) - Pj(E{j{Q))))\\ 2 be the statistic 
given by (3.10). We can derive now the following theorem as a special case 
of Theorem 3.1(a). 

Theorem 4.1. Assume j is the Veronese-Whitney embedding of M.P N ~ 1 
and {[X r ],\\X r \\ = l,r = 1, ...,n} is a random sample from a j-nonfocal 
probability measure Q on MP^ -1 that has a nondegenerate j-extrinsic vari- 
ance. Then T([v]) is given by 

(4.4) r(M) = ^*[K)a=i,...,iv-i]G(i,x)- 1 [( i , a ) 0=li ..., iV _ 1 ]* i /, 

and, asymptotically, T([u]) has a X%-i distribution. 

Proof. Since j is an isometric embedding and the tangent space Tj^jRP^" 1 
has the orthobasis u\,..., fjv-ij if we select the first elements of the adapted 
moving frame in Theorem 3.1 to be e a (Pj(vE,j)) = (d[v N ]j)(va), then the ath 
tangential component of Pj(j(X)) — Pj(v) w.r.t. this basis of T Pj ^ E (j(Q)))j(^-P N ~ 1 ) 
equals up to a sign the ath component of m — v^j w.r.t. the orthobasis 
ui, ... , vn-i in T^jMi?^ -1 , namely ^m. The result follows now from The- 
orem 3.1(a). □ 

Remark 4.1. If we apply Theorem 3.1(b) to the embedding j, we obtain 
a similar theorem due to Fisher, Hall, Jing and Wood [19], where T([^]) is 
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replaced by T([m]). Similar asymptotic results can be obtained for the large 
sample distribution of Procrustes means of planar shapes, as we discuss be- 
low. Recall that the planar shape space M = °f an ordered set of k points 
in C at least two of which are distinct can be identified in different ways 
with the complex projective space CP k ~ 2 (see, e.g., [8, 31]). Here we regard 
CP k ~ 2 as a set of equivalence classes CP k ~ 2 = S 2k ~~ 3 / S 1 where S 2k ~ 3 is the 
space of complex vectors in C k ~ l of norm 1, and the equivalence relation on 
g2k-3 - ls mu ltiplication with scalars in S l (complex numbers of modu- 
lus 1). A complex vector z = (z 1 , z 2 , . . . , z k_1 ) of norm 1 corresponding to a 
given configuration of k landmarks, with the identification described in [8], 
can be displayed in the Euclidean plane (complex line) with the superscripts 
as labels. If, in addition, r is the largest superscript such that z r ^ 0, then 
we may assume that z r > 0. Using this representative of the projective point 
[z] we obtain a unique graphical representation of [z], which will be called 
the spherical representation. 

The Veronese-Whitney (or simply Veronese) map is the embedding of 
CP k ~ 2 in the space of Hermitian matrices S(k — 1,C) given in this case by 
j ([z]) = zz* , where, if z is considered as a column vector, z* is the adjoint of 
z, that is, the conjugate of the transpose of z. The Euclidean distance in the 
space of Hermitian matrices S(k-l,C) is a%{A, B) = Tt((A - B){A - B)*) = 
Tr((A-B) 2 ). 

Kendall [31] has shown that the Riemannian metric induced on j(CP k ~ 2 ) 
by do is a metric of constant holomorphic curvature. The associated Rie- 
mannian distance is known as the Kendall distance and the full group of 
isometries on CP k ~ 2 with the Kendall distance is isomorphic to the special 
unitary group SU(/c — 1) of all (k — 1) x (k — 1) complex matrices A with 
A*A = I and det(A) = 1. 

A random variable X = [Z], \\Z\\ = 1, valued in CP k ~ 2 is j-nonfocal if the 
highest eigenvalue of E[ZZ*] is simple, and then the extrinsic mean of X is 
fJ'j E = Mi where v £ C fe_1 , ||zv|| = 1, is an eigenvector corresponding to this 
eigenvalue (see [8]). The extrinsic sample mean [z]j E °f a ran dom sample 
[z r ] = [(z^., . . . , z k_1 )], \\z r || = 1, r = 1, . . . ,n, from such a nonfocal distribution 
exists with probability converging to 1 as n — > oo, and is the same as that 
given by 

(4-5) ¥]j,e = H, 

where m is a highest unit eigenvector of 

n 

(4.6) K:=rT 1 ^z r z;. 

r=l 

This means that [2] ■ E is the full Procrustes estimate for parametric fami- 
lies such as Dryden-Mardia distributions or complex Bingham distributions 
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for planar shapes [35, 36]. For this reason, fij t E = [m] will be called the 
Procrustes mean of Q. 

Proposition 4.3. Assume X r = [Z r ],\\Z r \\ = l,r = 1, ... ,n, is a ran- 
dom sample from a j-nonfocal probability measure Q with a nondegenerate 
j -extrinsic covariance matrix on CP k ~ 2 . Then the j -extrinsic sample co- 
variance matrix G(j,X) as a complex matrix has the entries 

G(j, X) ab = n~ l (n k _i - r] a y 1 {ri k ^ 1 - n b y l 
( 4 ' 7 ) x ^2(m a ■ Z r ){m b ■ Z r )*|m fc _i • Z r \ 2 . 

r=l 

The proof is similar to that given for Proposition 4.2 and is based on the 
equivariance of the Veronese- Whitney map j w.r.t. the actions of SU (k — 1) 
on CP k ~ 2 and on the set S+(k — l, C) of nonnegative semidefinite self-adjoint 
(k — 1) by (k — 1) complex matrices (see [8]). Without loss of generality 
we may assume that K in (4.6) is given by K = diag{77 a } 0= i i ... i jt_i and 
the largest eigenvalue of K is a simple root of the characteristic polyno- 
mial over C, with m k —i = ek-\ as a corresponding complex eigenvector of 
norm 1. The eigenvectors over R corresponding to the smaller eigenvalues 
are given by m a = e a , m' a = ie a , a = 1, . . . , k — 2, and yield an orthobasis for 
T[m k _ 1 }j( < CP k ~ 2 ). For any z G S 2 ^ 1 which is orthogonal to mk-i in C^'" 1 
w.r.t. the real scalar product, we define the path 72(f) = [cosfmfc_i + sinfz]. 
Then T Pj ^j(CP k ~ 2 ) is generated by the vectors tangent to such paths 
72(f) at t = 0. Such a vector, as a matrix in S(k — 1,C), has the form 
zvn* k _ x + mk-iz* . In particular, since the eigenvectors of K are orthogo- 
nal w.r.t. the complex scalar product, one may take z = m a , a = 1, . . . , k — 2, 
orz = im a , a = 1, . . . ,k — 2, and thus get an orthobasis in T P .^j{M). When 
we norm these vectors to have unit lengths we obtain the orthonormal frame 

e a (Pj(K)) = d [rnkl] j(m a ) = 2~ 1/2 (m a m* k _ 1 + m fc _ira*), 

e' a ( P j( K )) = d[fn fc _i]i(»n»o) = i2~ 1/2 {m a m* k _ 1 - m k _ im * a ). 

Since the map j is equivariant we may assume that K is diagonal. In this 
case m a = e a , e a {Pj{K)) = 2" 1 / 2 ^- 1 and e' a (Pj(K)) = 2~ 1 / 2 F k ~\ where 
E\ has all entries zero except for those in the positions (a, b) and (6, a) that 
are equal to 1, and F% is a matrix with all entries zero except for those 
in the positions (a, b) and (6, a) that are equal to i, respectively —i. Just 
as in the real straightforward computation shows that dxPjiEa) = 

d K Pj(F^) = 0,Va < b < k - 1, and 

dKPjiEt 1 ) = (%_! - TlJ^eaiPjiK)), 
dKP^Ft 1 ) = (Vk-l ~ VaT^aiPjiK)). 
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We evaluate the extrinsic sample covariance matrix G(j,X) given in (3.8) 
using the real scalar product in S(k — 1,C), namely, U -V = ReTr([/V*). 
Note that 

dKPjiE^ 1 ) ■ e a (Pj(K)) = (%_! - jfo)" 1 ^, 
d K P j (Et 1 )-e' a (P j (K)) = 

and 

d K P J (F£- 1 )-e a (P J (K))=0. 

Thus we may regard G(j,X) as a complex matrix noting that in this case 
we get 

G{j,X) ab = rT 1 ^-! - l q a )~ 1 ij]k-x - r/ 6 ) _1 
^ 4 ' 8 ' ) x 5Z(e a • Z r )(e b ■ Z r )*\e k -i ■ Z r \ 2 , 

r=l 

thus proving (4.7) when K is diagonal. The general case follows by equiv- 
ariance. We consider now the statistic 

T(W E , f , E )=n\\G(j,XyVHMPj6W)-P 3 (VE))f 

given in Theorem 3.1 in the present context of random variables valued in 
complex projective spaces to get: 

Theorem 4.2. Let X r = [Z r ], \\Z r \\ = 1, r = 1, . . . ,n, be a random sam- 
ple from a Veronese-nonfocal probability measure Q on <CP k ~ 2 . Then the 
quantity (3.10) is given by 

(4.9) T([m], M) = n[(m ■ v a )a=i,...,k-2]G{j,X)- 1 [(m • v a ) a =i,...,k-2]* 
and asymptotically T([m], [u]) has a X2&-4 distribution. 

Proof. The tangent space T^ Uk _^CP k ~ 2 has the orthobasis u±, . . . , v k -2, v \- 
Note that since j is an isometric embedding, we may select the first el- 
ements of the adapted moving frame in Corollary 3.1 to be e a (Pj(/i)) = 
(d^ijiX^a), followed by e* a {Pj(^)) = (d^^j)^*). Then the ath tangen- 
tial component of Pj(j(X)) — Pj(fi) w.r.t. this basis of T P .^j(CP h ~ 2 ) equals 
up to a sign the component of m — v k —\ w.r.t. the orthobasis vi, . . . , v k -2 m 
T^ fe i ]CP fc_2 , which is z/*m; and the a*th tangential components are given 
by v^m, and together (in complex multiplication) they yield the complex 
vector [(m • fo)o=i,...,fc-2]- The claim follows from this and from (4.3), as a 
particular case of Corollary 3.1. □ 

We may derive from this the following large sample confidence regions. 
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Corollary 4.1. Assume X r = [Z r ], \\Z r \\ = 1, r = 1, . . . ,n, is a random 
sample from a j-nonfocal probability measure Q on CP k ~ 2 . An asymptotic 
(1 — a) -confidence region for fJ. J E (Q) = [v] is given by R a (K) = {[u] : T([m], [v]) < 
X2k-4a}' where T(\m\,[v\) is given in (4.9). If Q has a nonzero absolutely 
continuous component w.r.t. the volume measure on CP k ~ 2 , then the cover- 
age error of R a (X.) is of order 0(n _1 ). 

For small samples the coverage error could be quite large, and a bootstrap 
analogue of Theorem 4.2 is preferable. 

Theorem 4.3. Let j be the Veronese embedding of CP k ~ 2 , and let 
X r = [Z r ], \\Z r \\ = 1, r = l,...,n, be a random sample from a j-nonfocal 
distribution Q on CP k ~ 2 having a nonzero absolutely continuous component 
w.r.t. the volume measure on CP k ~ 2 . Assume in addition that the restric- 
tion of the covariance matrix of j(Q) to T^j(CP k ~ 2 ) is nondegenerate. Let 
^e(Q) = M be the extrinsic mean of Q. For a resample {^*}r-=i,...,n from 
the sample consider the matrix K* :=n~ x ^2iZ*Z* . Let (^)o=i,...,fe-i be the 
eigenvalues of K* in their increasing order, and let (m*) 0= i ... be the 
corresponding unit complex eigenvectors. Let G*(j,X)* be the matrix ob- 
tained from G(j,X) by substituting all the entries with * -entries. Then the 
bootstrap distribution function of 

T([m}*, [m]) := ^(m^ • m* a ) a = lr ..^ 2 ]G*({j,X)*y 1 [(rn k - 1 ■ m*) a= i r .. )fe _ 2 ]* 

approximates the true distribution function of T([m], [z/]) given in Theo- 
rem 4.2 with an error of order O p (n~ 2 ). 

Remark 4.2. For distributions that are reasonably concentrated one 
may determine a nonpivotal bootstrap confidence region using Corollary 
3.1(a). The chart used here features affine coordinates in CP k ~ 2 . Recall that 
the complex space C k ~ 2 can be embedded in CP k ~ 2 , preserving collinearity. 
Such a standard affine embedding, missing only a hyperplane at infinity, is 
(z 1 , . . . , z k ~ 2 ) — > [z 1 : ■ ■ ■ : z k ~ l : 1]. This leads to the notion of affine coordi- 
nates of a point 



p=[z L :---:z m :z k - 1 ], z^^O, 



to be defined as 



(w\w 2 ,...,w k - 2 ) = ( 



z 1 z k 2 



z k 1 ' ' z k 



To simplify the notation the simultaneous confidence intervals used in the 
next section can be expressed in terms of simultaneous complex confidence 
intervals. If z = x + iy,w = u + iv,x<u,y<v, then we define the complex 
interval (z, w) = {c = a + ib\a E (x, u),b £ (y, v)}. 
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5. Applications. In this last section we consider three applications. 

Application 1. Here we consider the data set of n = 50 South magnetic 
pole positions (latitudes and longitudes), determined from a paleomagnetic 
study of New Caledonian laterities ([20], page 278). As an example of appli- 
cation of Section 2, we give a large sample confidence region for the mean 
location of the South pole based on this data. The sample points to a non- 
symmetric distribution on S 2 ; the extrinsic sample mean and the intrinsic 
sample mean are given by 

X E = (0.0105208, 0.199101, 0.979922)* 

and, using X E as the initial input of the necessary minimization for con- 
structing Xj, 

X I=p = (0.004392, 0.183800, 0.982954)* . 

From Examples 2.1 and 2.2, select the orthobasis ei(p),e2(p) given in (2.3) 
and the logarithmic coordinates u l ,u 2 w.r.t. this basis in T p S 2 defined in 
(2.4). Then compute the matrix A given in (2.22), to get, using Corollary 
2.2, the following 95% asymptotic confidence region for fj,j: 

U = {Exp p (u 1 ei(p) + u 2 e 2 {p))\ 

16.6786(n 1 ) 2 - 2.9806uV + lO^lSO^ 1 ) 2 < 5.99146}. 

Note that Fisher, Lewis and Embleton ([20], page 112) estimate another 
location parameter, the spherical median. The spherical median here refers 
to the minimizer of the expected geodesic (or, arc) distance to a given point 
on the sphere. For this paleomagnetism data, their sample median is at 
78.9°, 98.4°, while the extrinsic sample mean is 78.5°, 89.4° and the intrinsic 
sample mean is 79.4°, 88.6°. These estimates differ substantially from the 
current position of the South magnetic pole, a difference accounted for by 
the phenomenon of migration of the Earth's magnetic poles. 

Application 2. As an application of Section 4, we give a nonpivotal 
bootstrap confidence region for the mean shape of a group of eight landmarks 
on the skulls of eight-year-old North American children. The sample used 
is the University School data ([10], pages 400-405). The data set represents 
coordinates of anatomical landmarks, whose names and position on the skull 
are given in [10]. The data are displayed in Figure 1. (The presentation of 
raw data is similar to other known shape data displays such as in [15], page 
46.) The shape variable (in our case, shape of the eight landmarks on the 
upper mid face) is valued in a planar shape space CP 6 (real dimension = 12). 
A spherical representation of a shape in this case consists of seven marked 
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University School Data - Boys 
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Fig. 1. 

points; in Figure 2 we display a spherical representation of this data set. A 
representative for the extrinsic sample mean (spherical representation) is 

(-0.67151 + 0.66823i, 0.76939 + 1.05712i, -1.03159 - 0.15998i, 

-0.57776 - 0.87257^,0.77871 - 1.36178?, 

-0.17489 + 0.82106i, 1.00000 + O.OOOOOi). 

We derived the nonpivotal bootstrap distribution using a simple program in 
S-Plus4.5, that we ran for 500 resamples. A spherical representation of the 
bootstrap distribution of the extrinsic sample means is displayed in Figure 
3. Here we added a representative for the last landmark (the opposite of the 
sum of the other landmarks since data is centered at 0). 
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University School Data - Spherical Representation 




Note that the bootstrap distribution of the extrinsic sample mean is very 
concentrated at each landmark location. This is in agreement with the the- 
ory, that predicts in our case a spread of about six times smaller than the 
spread of the population. It is also an indication of the usefulness of the 
spherical coordinates. We determined a confidence region for the extrinsic 
mean using the six 95% simultaneous bootstrap complex intervals for the 
affine coordinates, as described in Remark 4.2, and found the following com- 
plex intervals: 



R. BHATTACHARYA AND V. PATRANGENARU 



Bootstrap Distribution of 500 Extrinsic 
Sample Mean Configurations 



s e 



f 

I I I 1 I 

-0.4 -0.2 0.0 0.2 0.4 

X 

Fig. 3. 

(-0.677268 + 0.666060?, -0.671425 + 0.672409i), 

(0.767249 + 1.051660i,0.775592 + 1.058960i), 
(-1.036100 - 0.161467z, -1.029420 - 0.154403i), 
(-0.578941 - 0.875168i, -0.574923 - 0.871553i), 
(0.777688 - 1.366880^,0.782354 - 1.358390i), 
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for wq: 



(-0.177261 + 0.820107i, -0.173465 + 0.824027i). 



Application 3. This example is relevant in glaucoma detection. Al- 
though it is known that increased intraocular pressure (IOP) may cause a 
shape change in the eye cup, which is identified with glaucoma, it does not 
always lead to this shape change. The data analysis presented shows that 
the device used for measuring the topography of the back of the eye, as 
reported in [11], is effective in detecting shape change. 

We give a nonpivotal bootstrap confidence region for the mean shape 
change of the eye cup due to IOP. Glaucoma is an eye disorder caused by 
IOP that is very high. Due to the increased IOP, as the soft spot where 
the optic nerve enters the eye is pushed backwards, eventually the optic 
nerve fibers that spread out over the retina to connect to photoreceptors 
and other retinal neurons can be compressed and damaged. An important 
diagnostic tool is the ability to detect, in images of the optic nerve head 
(ONH), increased depth (cupping) of the ONH structures. Two real data- 
processed images of the ONH cup surface before and after the IOP was 
increased are shown in Figure 4. 

The laser image files are, however, huge-dimensional vectors, and their 
sizes usually differ. Even if we would restrict the study to a fixed size, there 
is no direct relationship between the eye cup pictured and the coordinates 
at a given pixel. A useful data reduction process consists in registration 
of a number of anatomical landmarks that were identified in each of these 
images. Assume the position vectors of these landmarks are X\, . . . , Xk, k > 
4. Two configurations of landmarks have the same shape if they can be 
superimposed after a translation, a rotation and a scaling. The shape of the 
configuration x= (xi, . . . ,x^) is labelled o(x) and the space of shapes 



Fig. 4. Change in the ONH topography from normal (left) to glaucomatous (right). 
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of configurations of k points in R m at least two of which are distinct is the 
shape space introduced by Kendall [31]. 

We come back to the shape of an ONH cup. This ONH region resembles a 
"cup" of an ellipsoid and its border has a shape of an ellipse. In this example 
four landmarks are used. The first three landmarks, denoted by S, T and 
N, are chosen to be the "top, left and right" points on this ellipse, that is 
(when referring to the left eye), Superior, Templar and Nose papilla. The 
last landmark V that we call vertex is the point with the largest "depth" 
inside the ellipse area that determines the border of the ONH. Therefore, in 
this example the data analysis is on the shape space of tetrads S3, which 
is topologically a five-dimensional sphere (see [33], page 38); however, the 
identification with a sphere is nonstandard. On the other hand, it is known 
that if a probability distribution on has small support outside a set of 
singular points, the use of any distance that is compatible with the orbifold 
topology considered is appropriate in data analysis ([15], page 65) since 
the data can be linearized. Our choice of the Riemannian metric (5.3) is 
motivated by considerations of applicability of Theorems 2.2 and 2.3 and 
computational feasibility. Dryden and Mardia ([15], pages 78-80) have in- 
troduced the following five coordinates defined on the generic subset of £3 
of shapes of a nondegenerate tetrad that they called Bookstein coordinates: 

v 1 = (w 12 w 13 + w 22 w 23 + ^32^33)/^, 

V 2 = {(W 12 W 2 3 - W22Wr3,) 2 



V 



V 



A 



,3 



+ (^12^33 - W 32 W 13 ) 2 + (^22^33 ~ ^ 23^32) 2 ) 1/2 / a ; 
(tVuWu + W 22 W 2 4 + W 32 Wu)/a, 

{ab l/2 )~ l (w\ 2 (w 23 w 2 4 + w 33 w u ) + wl 2 (w 13 w 14: + W33W34) 



(5.1) 



+ ^32(^13^14 + ^23^24) - Wl 2 Wl 3 (w 22 W 2 4 + W 32 W 3 4) 
~ W 22 W 3 2(W23W34 + W 33 W 2 4:) 

~ W 12 Wu(w22W2 3 + W 32 W 33 )) 



V 



,5 



(^12^23^34 - ^12^33^24 ~ ^13^22^34 
+ ^13^32^24 + ^22^33^14 - ^32^23^14)/ (2o6) 



where 



a = 2{w\ 2 + wl 2 + wl 2 ), 
(5.2) b = w\ 2 w 23 + w\ 2 w 33 - 2wi 2 wi 3 w 22 w 23 + w\ 3 w 22 + u> 2 3 u>§ 2 

- 2lUi2U;i3W32W33 + ^33^22 + ^23^32 ~ 2^22^32^23^33 



and 



Wri =X\ - (x r 1 +X 2 )/2, 



r = 2,3,4. 
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These coordinates carry useful geometric information on the shape of the 
4-ad; v 1 and t> 3 give us information on the appearance with respect to the bi- 
sector plane of [XiA2],f 2 and v 4 give some information about the "flatness" 
of this 4-ad and v 5 measures the height of the 4-ad (Xi, X2, A3, X4) relative 
to the distance \\Xi — X^W- Assume U is the set of shapes o(X) such that 
(X\, X2, A3, Xi) is an affine frame in M 3 , and <j>:U -> R 3k ~ 7 is the map that 
associates to o(X) its Bookstein coordinates. U is an open dense set in S3, 
with the induced topology. In the particular case k = 4, S3 is topologically 
a five-dimensional sphere and, from a classical result of Smale [48], £3 has 
a differentiable structure diffeomorphic with the sphere S 5 . Moreover, if L 
is a compact subset of U, there are a finite open covering U\ = U, . . . ,Ut of 
S3 and a partition of unity 4>i, . . . , <pt, such that 4>i(o(X)) = 1, Vo(A) £ L. 
We will use the following Riemannian metric on S|: let (yi, ■ • ■ , 2/5) be the 

Bookstein coordinates of a shape in U\ and let g\ = dy\ H Y dy\ be a flat 

Riemannian metric on U\, and for each j = 2, . .. ,t we consider any fixed 
Riemannian metric gj on Uj. Let g be the Riemannian metric given by 

t 

(5.3) 9 = ^24>j9j- 

i=y 

The space (S|,d g ) is complete and is flat in a neighborhood of L. In this 
example the two distributions of shapes of tetrads before and after increase 
in IOP are close. Hence L, which contains supports of both distributions, 
consists of shapes of nondegenerate tetrads only. 

Computations for the glaucoma data yield the following results. The p- 
value of the test for equality of the intrinsic means was found to be 0.058, 
based on the bootstrap distribution of the chi square-like statistic discussed 
in Remark 2.6. The number of bootstrap resamples for this study was 3000. 
The chi square-like density histogram is displayed in Figure 5. A matrix plot 
for the components of the nonpivotal bootstrap distribution of the sample 
mean differences 7* in Remark 2.6 for this application is displayed in Figure 
6. The nonpivotal bootstrap 95% confidence intervals for the mean differ- 
ences 7j, j = 1, . . . , 5, components of 7 in Remark 2.6 associated with the 
Bookstein coordinates Vj, j = 1, . . . , 5, are: (—0.0377073, —0.0058545) for 71, 
(0.0014153,0.0119214) for 72 , (-0.0303489,0.0004710) for 73, (0.0031686,0.0205206) 
for 74, (-0.0101761,0.0496181) for 75. Note that the individual tests for 
difference are significant at the 5% level for the first, second and fourth 
coordinates. However, using the Bonferroni inequality, combining tests for 
five different shape coordinates each at 5% level leads to a much higher 
estimated level of significance for the overall shape change. 

APPENDIX 

The data set in Application 3 consists of a library of scanning confocal 
laser tomography (SCLT) images of the complicated ONH topography [11]. 
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Fig. 5. x 2 -tike bootstrap distribution for equality of intrinsic mean shapes from glaucoma 
data. 
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Bookstein coordinates due to increased IOP. 



INTRINSIC AND EXTRINSIC MEANS— II 



33 



Those images are the so-called range images. A range image is, loosely speak- 
ing, like a digital camera image, except that each pixel stores a depth rather 
than a color level. It can also be seen as a set of points in three dimensions. 
The range data acquired by 3D digitizers such as optical scanners commonly 
consist of depths sampled on a regular grid. In the mathematical sense, a 
range image is a 2D array of real numbers which represent those depths. 
All of the files (observations) are produced by a combination of modules in 
C++ and SAS that take the raw image output and process it. The 256 x 256 
arrays of height values are the products of this software. Another byproduct 
is a file which we will refer to as the "abxy" file. This file contain the fol- 
lowing information: subject names (denoted by: lc, Id, le, If, lg, li, lj, lk, 
11, In, lo, lp), observation points that distinguish the normal and treated 
eyes and the 10° or 15° fields of view for the imaging. The observation point 
"03" denotes a 10° view of the experimental glaucoma eye, "04" denotes a 
15° view of the experimental glaucoma eye, "11" and "12" denote, respec- 
tively, the 10° and the 15° view of the normal eye. The two-dimensional 
coordinates of the center (a, b) of the ellipses that bound the ONH region, 
as well as the sizes of the small and the large axes of the ellipses (x,y), are 
stored in the so-called "abxy" file. To find out more about the LSU study 
and the image acquisition, see [11]. File names (each file is one observa- 
tion) were constructed from the information in the "abxy" file. The list of 
all the observations is then used as an input for the program (created by 
G. Derado in C++) which determines the three-dimensional coordinates of 
the landmarks for each observation considered in our analysis, as well as for 
determining the fifth Bookstein coordinate for each observation. Each image 
consists of a 256 x 256 array of elevation values which represent the "depth" 
of the ONH. By the "depth" we mean the distance from an imaginary plane, 
located approximately at the base of the ONH cup, to the "back of the ONH 
cup." 

To reduce the dimensionality of the shape space to 5, out of five landmarks 
T, S, N, I, V recorded, only four landmarks {X\ = T, X^ = S, A3 = N, 
X4 = V) were considered. 

The original data were collected in experimental observations on Rhesus 
monkeys, and after treatment a healthy eye slowly returns to its original 
shape. For the purpose of IOP increment detection, in this paper only the 
first set of after-treatment observations of the treated eye is considered. 
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